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ABSTRACT OF THE DISCLOSURE 

A method and apparatus for the automatic extraction of 
linear features from digital imagery is disclosed herein. The proposed 

5 methodology derives from the Radon transform, yet has none of its basic 
disadvantages. Indeed, it is not limited solely to extracting straight lines, 
there are no false peaks or virtual maxima (peaks that do not represent 
the real line on the input image) in the transformed plane and the process 
does not have difficulty detecting line segments which are significantly 

1 0 shorter than image dimensions. 
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TITLE OF THE INVENTION 

Automatic Extraction of Linear Features from Digital Imagery 

5 

FIELD OF THE INVENTION 

The present invention relates to digital imagery. More 
10 specifically, the present invention is concerned with the automatic 
extraction of linear features from digital imagery. 

BACKGROUND OF THE INVENTION 

15 

With the recent, rapid development of computer tools, 
digital image processing has become an efficient, economic, and 
convenient means for obtaining qualitative and quantitative information in 
different fields such as remote sensing, cartography, robotics, and 

20 materials. Indeed, digital image processing makes it possible to 
qualitatively describe images from various sources in terms of pattern 
recognition to identify and isolate contained objects. One of the major 
subclasses of contained objects often present in digital images is the 
linear features. Automatic detection of linear features from digital images 

25 plays an important role in pattern recognition and digital image processing. 

A variety of techniques and many algorithms have 
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emerged to automatically extract linear features from digital images. 
These techniques can be classified into two main categories: local 
methods (which are, based on local operators such a mobile kernel) and 
global methods (which focus on mathematical transformations, such as 
5 Hough transforms) [1]. 

Local methods for automatically extracting linear 
features exploit local variations of pixel intensity in a small neighborhood 
by calculating of the gradients in small, limited-size windows in the image, 

10 e.g., 3 x 3 or 5 x 5 [2][3][4][5][6][7][8][9]. A number of researchers have 
examined mathematical morphology as a means of extracting linear 
features [10][11][12][13][14]. The problems pertaining to this technique 
arise from the number of human decisions required to reconnect and 
rebuild line segments, which increases processing time. Multi-dimensional 

1 5 line detection is the other technique for detecting for linear features that 
collects different spectral information for the same scene and may 
highlight different parts of lines [15]. The first stage of this method for 
obtaining the combined images requires several transformations of 
multiple original bands. Human intervention is needed to select the best- 

20 combined image. Another approach to linear-feature extraction involves 
knowledge-based systems, which need more information than a simple 
digital image for line extraction [16][17][18]. 

These local methods generally remain inefficient 
25 because they fail to have a global view of the linear features in a digital 
image. One problem common to all of these methods is that the resulting 
extracted line images contain a fair amount of noise, while the detected 
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lines are incomplete and geometrically shifted. These difficulties are 
magnified by intersections and linear features that display some curvature 
[19][20][1]. In addition, these methods turn in exceptionally long 
processing times when extracting features from large images [16]. 

5 

The Radon transform and its derivative, the Hough 
transform, are the most frequently used approaches as global methods for 
detecting linear features [21][22][23][24]. In principle, a straight line from 
the input image is transformed into a digital peak (a light or dark pixel, 
1 0 compared to its neighborhood) in the transformed plane. In this case, it is 
easier to detect a peak in the transformed plane than straight line 
detection in the input image. There are three basic limitations for these 
methods that sorely restrict their applications and their utility in practice. 

1 5 Firstly, the Hough transform-based methods are limited 

solely to extracting straight lines [21]. Therefore, linear features that span 
the entire image but display some curvature may not produce suitable 
peaks or troughs in the transform plane. This restriction is linked directly 
to the basic definition of the method. 

20 

Secondly, there are false peaks or virtual maxima (peaks 
that do not represent the real line on the input image) in the transformed 
plane. These false peaks considerably lower the quality of the results by 
increasing the error of commission. When the line density in the input 
25 image is high, eliminating the false peaks from transformed plane 
seriously limits these methods. 
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The last methodological inconvenience of these 
approaches is the ignorance of the nature of the detected lines. Since 
intensity integration in the transformation process is performed over the 
entire length of the image, the process can have difficulty detecting line 
5 segments which are significantly shorter than image dimensions 
[21][25][26]. Neither can it provide information about the positions of the 
endpoints of these shorter line segments or line length. 

10 OBJECTS OF THE INVENTION 

An object of the present invention is therefore to provide 
an improved approach for the extraction of linear features. 

15 Other objects, advantages and features of the present 

invention will become more apparent upon reading of the following non- 
restrictive description of preferred embodiments thereof, given by way of 
example only with reference to the accompanying drawings. 

20 

BRIEF DESCRIPTION OF THE DRAWINGS 

In the appended drawings: 

25 Figures 1a, 1b and 1c are a flowchart and a legend of 

the proposed technology; 
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Figure 2 is a binary image simulated to showcase the 
method's ability to extract all lines without consideration for the slope 
line; extract the lines having a shorter length of the image dimension; 
and locate extracted lines accurately; 

5 

Figure 3 is a binary image 2 obtained by applying the 
implemented algorithm to the binary image of Figure 2; 

Figure 4 is a binary image simulated in order to 
10 demonstrate the method's performance in extracting curved as well as 
straight lines; 

Figure 5 is a binary image produced by applying the 
implemented algorithm to the binary image of Figure 4; 

15 

Figure 6 is a binary image of a part of a digitized road 
map of the city of Sherbrooke; and 

Figure 7 is a binary image 6 produced by applying the 
20 implemented algorithm to the binary image of Figure 6. 



DESCRIPTION OF THE PREFERRED EMBODIMENT 



25 



The proposed methodology derives from the Radon 
transform, yet has none of its basic disadvantages already mentioned 
hereinabove. It also builds on a solid mathematical base which 
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demonstrates the conceptual feasibility of the new method 
mathematically. An algorithm has also been developed and implemented 
in order to achieve a concrete realization of the innovative method in a 
real context. 



as general procedures for a algorithmic development of our proposed 
method will be presented in detailed. 

1 0 Mathematical Development 

Continuous function 

Let g(x,y) be a continuous signal of the continuous 
15 variables x and y and let £ denote a ri-dimensional parameter vector 
defined as: 



5 



In following subsections, the mathematical base as well 



(D 



20 



where £ spans the parameter domain. 



For a two-dimensional continuous function, the 



parameter vector can be defined as: 



25 



Tl = 3 



(2) 



or: 
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S = (a,P,Y) (3) 

hence: 

g(x 5 y) = <t>(x 5 y;^) = y-ax 2 -Px-y (4) 

5 The Radon transform g(£) of function g(x 9 y) \s defined 

as: 

8® = £ £g(x,y)5(<Kx,y;S)dxdy (5) 

10 where 8(.) denotes the Dirac delta function. 

Using the definition (5) of Radon transform, curves 
expressed by the next parameter form can be detected: 

15 <Kx,y;$) = 0 (6) 

Substituting the relation (4) into equation (5), we find: 
g(a,p,y)= £ £g(x,y)8(y-ax 2 -px-y)dxdy (7) 

20 

where y = <t>(x;£) represent the transformation curve. 

Using proprieties of delta function, equation (7) 

becomes: 



25 
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g(a,p,y)= ££g(x,y)8(y-ax 2 - Px - y)dxdy = £g(x,ax 2 + px + y)dx 



(8) 



At this stage, the transformation curve is defined as a 
polynomial function of second degree using the delta function as below: 



g(x,y) = S(y-<Kx£*)) 



(9) 



where the parameter vector of this curve is: 
ft'-(a\P\T) 

Using equations (4), (9) and (10), we find: 
g(x,y) = 5(y-(a*x 2 +p*x+y*)) 



(10) 



(11) 



According to definition of the Radon transform: 
g(©=£8tt(x;9-*(x;S*))dx 

= L±T 



8(x 




a*(x;y 


awx;5*) 


ax 


3x 



a*(x,;0 


aKx,**) 


ax 


ax 



(12) 
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Substitution of parameter vectors (3) and (10) into the 
transform definition (12) gives: 

g(ct,p, y) = £6((ax 2 + Px + y) - (a * x 2 + p * x + y*)) dx 



6(x-x a ) 



d(ax 2 +px + y) g(a*x 2 + p*x + y*] 
fix dx 



,dx (13) 



10 



g(ax 2 + px + y) d(a*x 2 +p*x + y*) 
dx dx 



Finally, we find: 



g(a,p,y)=oo 



pour a=a* et P = p* et y = y* 



0 -< g(a,p,y) oo pour autres cas 
15 (14) 



Equation (14) shows that the Radon transform of a 
polynomial function of second degree will give an infinite value (a peak) 
in the parameter domain, when parameters of this function are precisely 
20 determined. In addition, this equation confirm that the parameter domain 
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will contain some other non-zero values because of parameters values 
which have some similarity with real parameter values of the polynomial 
function in question. 

5 In the light of the mathematical demonstration presented 

above, tree important issues are concluded: 

i) The Radon transform of a two-dimensional plane containing a 
polynomial function of second degree is a Euclidean tree- 
dimensional space; 

1 0 ii) this polynomial function is mapped into a unique point in the tree- 
dimensional space; 
iii) position of this point in the tree-dimensional space determines the 
polynomial function parameters. 

15 The following subsection provides a discrete 

approximation to the Radon transform for a polynomial function of second 
degree. This discrete approximation is the key of the algorithmic 
development of the polynomial function in question. 

20 Discrete Function 

Discrete Radon transform of a discrete function 
4(m;l,k,h) is defined as: 

M-l M-l 

25 g(l,k,h) = £g(m,<Km;l,k,h))= £g(m,n(m;l,k,h)) (15) 

m=0 ro*0 

where : 
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<j>(m;l ) k,h)=n = 
and 



4>(x min +mAx;e(l,k,h)-y mia ) 
Ay 



(16) 



x =x .„+mAx 



m rain 



m = 0,l,-,M 
n = 0,1,— ,N 



(17) 



where M and N are the discrete plane dimensions and Ax, Ay are sampling 
intervals of this plane. 

For a polynomial function of second degree expressed 
as y - ex 2 - px - x = 0 , the Radon transform is: 

g(l,k,h) = Ax^gta,^ + Pk x m +t h ) (18) 

with: 



£l =e rain + IAe 

P lt =P,nin+ kAP 

l*h=*mm + ^ 



1 = 0,1, -,L 
k = 0,l,-".K 
h = 0,l,--,H 



(19) 



where L, K and H represent dimensions of Radon domain and Ae, Ap and 
At are sampling intervals of this domain. 



Using relation (19), the function n(m;l,k,h) of relation 

(16) becomes: 



CA 02313803 2000-07-11 



12 



n(m;l,k,h) = 



r S|Xj+P fc x.+T| 
L Ay 



h / roia 



(20) 



With substituting the relation (17) in the previous 



equation, we find: 



n(m;l,k,h) = 



,(x min + mAx) 2 + p k (x min + mAx) + x h - y min 1 
Ay J 



(21) 



For simplifying, the previous equation can be expressed 



as: 



10 



n* = am 2 + pm + y 



(22) 



where: 



15 



a = 



s,Ax 2 
Ay 



P^ ( 2 x mi ,g| + Pic)Ax 
Ay 

g|xL+PkX min +T h -y n 
Ay 



(23) 



According to the equation (19), once the Radon domain 
parameters has been initialized, using equations (21), (22) and (23) this 
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domain is determined for all of two-dimensional plane points (X,Y). 
Algorithmic development 

The concrete achievement of the proposed method is 
schematized in a flowchart presented in Figures 1a and 1b. This flowchart 
illustrates also all essential instructions for an algorithmic development of 
the proposed technology. 

1) Input image 

Represents the digital input image. 

2) Design the discrete parameter domain 

The parameter vector of the transformation is defined in this step as 
described in equations (1), (2) and (3) in subsection (3.1) (Mathematical 
Development). 

S = (a,p,Y) ( 3 ) 

3) Length of lines to be detected 

Input of the interval of line length (i.e. lines between 5 and 220 pixels). 
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4) Initialize parameter domain 

This step is for initializing the numerical values of discrete parameters 
8^ Ae, Pmin , Ap, Ai, L, K et H, as described in equation (19) and 
x m ,y n , Ax, Ay, as described in equation (17) (M and N correspond to the 
5 dimension of digital input image). 



K=x min +mAx 

|y„ = y^in + nty 



m = 0,l,",M 
n =0,1,-,N 



(17) 



10 



|£, =6 



+ IAe 

Pk = Pmin + 



1 = 0,1,-,L 
k = 0,l,-,K 
h = 0,l,-,H 



(19) 



5) Start of 1 st loop (for every position in the parameter domain) 

15 If s a loop for determining every possible combination of e„ p k et i„ as 
described in equation (19) by incrementing tree values of I, k and I (0 < I 
< L, 0 < k< K, 0 < h < H). 



20 6) Computation of line / curve deterministic parameters 
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In this step, the deterministic parameters a, p and y are computed by 
using equation (23). 



e,Ax 2 
a = — 



Ay 

o_ (2x min e,+p k )Ax 
Ay 

v = g l X min + Pk X min + T h ~ yroin 

Ay 



(23) 



7) Start of 2 nd loop (determination of corresponding pixel's co- 
ordinates) 

In this loop, for each value of m (0<m<M), a new parameter n* is 
10 computed by using equation (21). 



n(m;l,k,h) = [^ m ^mAx) 2 + p k ( X|nin + mAx) + x h 

L A y 



(21) 



15 8) Finding the nearest pixel to determined co-ordinates in the input 
image 

In this step, the determined co-ordinates (m, h*) of the previous step 
are used to find the nearest neighbour approximation pixel P(m, n) 
within the input image. 



20 
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9) Numerical value of pixel is zero 

Ifs a verification step for numerical value of determined pixel P(m, n) in 
step 8. 

5 10) Increment number of zero pixels 

If the numerical value of determined pixel P(m, n) in step 9 is zero, the 
number of zero pixels will be incremented. 

11) Save the pixel's co-ordinates in the 1 st layer of database 

1 0 The pixel's co-ordinates m and n are stored in the first layer of the 
database. 

12) Save the obtained value in the flayer of database 

The numerical value of the determined pixel P(m t n) in step 9 is stored 
1 5 in the 4 th layer of database. 

13) Increment number of non zero pixels 

If numerical value of the determined pixel P(m, n) in step 9 is not zero, 
the number of non zero pixels will be incremented. 

20 

14) Save the obtained value in the 3 rd layer of database 

The numerical value of the corresponding pixel is stored in the 4 th layer 
of the database. 

25 15) Save the pixel co-ordinates in the 2 nd layer of database 

The non-zero pixel's co-ordinates m and n are stored in the second 
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layer of the database. 

16) Accumulation of pixel numerical values 

In this step, the numerical values of pixels are accumulated. 

5 

17) Save the obtained value in the 5 th layer of database 

The value obtained in the step 16 is stored in the 5* layer of the database. 

18) End of the 2 nd loop 

1 0 This is the end of the 2 nd loop (loop for all values of m 0 < m <M). 

19) End of the 1 st loop: 

This is the end of the 1 st loop (loop for all values of I, k and h). 

1 5 20) Start of the 3 rt1 loop (for every cell in the 3 rt layer) 

It's a loop for recalling all the numerical values of non-zero pixels 
already stored in the 3" layer of the database. 

21) The value is inside the predefined range 

20 In this step, the determined numerical value of precedent step is 
verified to find if it's inside the predefined interval in step 3. 

22) Provided database 

It's the provided database. 

25 

23) Finding endpoints of the line 
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If the respond of step 21 is positive, stored co-ordinates of the 
corresponding pixels in the 3 rd layer of the database are recalled. Then, 
in the input image, the numerical values of these pixels will be verified 
to find endpoints of the corresponding line. 

5 

24) Restoration of determined line in output image 

In this step, by using the stored deterministic parameters of the line and 
the line endpoints that are found in the precedent step, the detected 
line will be restored in the output image. 

10 

25) Save ail the extracted line information in a new database 

All the information about extracted lines, such as endpoints, lengths, 
etc. are stored in a new database. 

15 26) End of the 3 rd loop 

This is the end of the 3 rd loop (loop for all cells in the 3 rd layer of the 
database). 

27) End of algorithm 

20 The end of algorithm. 

Features and Operation of the Invention 

We have developed a new global method that is capable 

25 of: 

- detecting and extracting the lines from a digital image with any 



CA 02313803 2000-07-11 



19 



5 



10 



curvature; 

- discarding the virtual maxima on the Hough or on the Radon 
transformed plane; 

- locating every line segment in the input image (endpoints and length). 

Demonstration and Discussion 

For a visual demonstration, the urban road network 
detection of the city of Sherbrooke from satellite imagery has been 
selected as an application. 



The binary image of Figure 2 was simulated to showcase 
the method's ability to: 

i) extract all lines without consideration for the slope line; 

ii) extract the lines having a shorter length of the image dimension; 
1 5 iii) locate extracted lines accurately. 



Figure 3 is a binary image obtained by applying the 
implemented algorithm to the image of Figure 2. The absolute similarity of 
Figure 3 with Figure 2 demonstrates clearly the accuracy of three abilities 
20 of the method as described above. 



The binary image of Figure 4 was simulated in order to 
demonstrate the method's performance in extracting curved as well as 
straight lines. The binary image illustrated in Figure 5 was produced by 
25 applying the implemented algorithm to the binary image of Figure 4. The 
last image demonstrates this performance. 
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The binary image of Figure 6 is part of a digitized road 
map of the city of Sherbrooke. Finally, the binary image of Figure 7 was 
produced by applying the implemented algorithm to the binary image of 
5 Figure 6. A comparison of the binary images of Figures 6 and 7 reveals 
that the restrained roads in the input image have been extracted in detail 
and with high geometric precision. 

Advantages and Disadvantages of the Invention versus Competitive 
10 Technologies 

To the best of our knowledge, there is no operational 
algorithm that automatically extracts line features and provides accurate, 
acceptable results either referred to in the literature or available on the 
15 market. The interactive and semi-automatic methods remain the most 
used. 

Of the local methods, LINDA (Linear-featured Network 
Detection and Analysis) system [2][7][27], based on the profile intensity 
20 analysis of the pixel line, is the most recent and the most representative 
method. With all of basic disadvantages of the local methods (cf. the first 
paragraph of page 5), the LINDA system is far from being operational with 
respect to systematically and automatically processing a large set of 
images. 

25 

As for the global method, the Hough transform is 
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probably the method most frequently mentioned and used by scientific 
researchers in various fields [9][21][25]. Due to its basic limitations as 
discussed above, the method is to a large degree limited to a specific type 
of application in which lines are straight, continuous, and of low density. 

When compared to techniques reported in the scientific 
literature or now commercially available, the proposed method is unique, 
since it can extract curved lines of any length with a high degree of 
accuracy. Our method identifies every line or line segment, labeling its 
length, its curvature, and endpoint coordinates in a database. Moreover, 
the database can be easily integrated into any geographic or spatial 
information system. Consequently, everything points to our method as 
having significant potential. It could lead to the development of robust, 
operational software to automatically extract line features from any kind of 
digital image from various sources. 

Examples of Fields of Applications 

The method is of broad interest and makes a general 
contribution to the field of pattern recognition and digital image processing, 
regardless of the nature or source of the digital image. It could be used in 
a vast number of fields and application, of which the following are a 
sampling. 

A) Remotely sensed imagery 

- Extracting and updating road maps. 
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- Planimetric cartography, geology, and hydrology. 

- Mapping power-lines, pipelines, and railways. 

- Tracking cycle trails, ORV trails, and the like. 

5 B) Medical imagery 

- Detecting of arteries, blood vessels, tumors, etc; 

C) Imagery from other source 

- Electrical circuits and the like. 

10 
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Although the present invention has been described 
hereinabove by way of preferred embodiments thereof, it can be modified, 
without departing from the spirit and nature of the subject invention as 
15 defined in the appended claims. 
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WHAT IS CLAIMED IS: 

1. A method for the automatic extraction of linear features 
from digital imagery substantially as disclosed herein. 

2. An apparatus for the automatic extraction of linear 
features from digital imagery substantially as disclosed herein. 
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